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ABSTRACT 

Linear  prediction  has  become  an  important  tool  for  stationary  time  series  analysis. 
The  all-pole  system  model  which  is  a  product  of  the  linear  predictive  approach  has  ap- 
plications in  numerous  engineering  problems.  This  thesis  develops  a  simple  method  for 
obtaining  the  two-dimensional  all-pole  system  model.  The  lattice  structures  that  can  be 
used  to  implement  the  prediction  error  and  synthesis  filters  are  also  shown  to  have  an 
analogous  two-dimensional  counterpart.  The  construction  of  these  filters  is  in  terms  of 
orthogonal  Szego  polynomials  which  can  be  used  to  solve  the  two-dimensional  block 
Toeplitz  normal  equations  in  a  recursive  manner.  This  recursion  not  only  leads  to  the 
two-dimensional  (2-D)  lattice  structure,  but  also  allows  for  expansion  of  the  filter  order 
without  resolving  the  normal  equations.  Several  examples  are  presented  using  the  two- 
dimensional  linear  prediction  results  for  spectral  estimation  and  signal  synthesis. 
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I.     INTRODUCTION 

A.     LINEAR  PREDICTION 

The  initial  goal  oflinear  prediction  is  to  estimate  the  present  value  of  a  random  se- 
quence based  on  a  linear  combination  of  all  the  past  values.  This  is  equivalent  to  de- 
termining a  linear,  shift  invariant,  causal  prediction  error  filter  which  whitens  the 
random  sequence.  There  has  been  a  great  deal  of  study  of  one-dimensional  linear  pre- 
diction analysis  and  the  usefulness  of  the  resulting  filters.  The  purpose  of  this  thesis  is 
to  extend  in  a  straightforward  manner  the  results  of  one-dimensional  linear  prediction 
to  two-dimensional  signals. 

There  are  several  efficient  methods  for  the  solution  of  the  normal  equations  that 
arise  in  one-dimensional  linear  prediction  problems.  The  most  familiar  of  these  methods 
is  the  Levinson  algorithm.  The  Levinson  algorithm  alleviates  the  need  to  perform  a 
matrix  inversion  in  the  solution  of  the  normal  equations  for  the  filter  coefficients.  The 
recursive  nature  of  the  algorithm  also  allows  the  prediction  error  filter  to  be  implemented 
in  a  lattice  structure.  These  two  features  will  be  of  primary  concern  in  the  solution  of 
the  two  dimensional  problem. 

There  are  two  general  approaches  to  linear  prediction.  The  first  method,  which  is 
usually  referred  to  as  the  spectral  factorization  method,  begins  with  the  knowledge  of  the 
power  spectral  density  of  the  random  process  and  attempts  to  find  the  model  parameters 
that  fit  the  spectrum.  This  works  well  in  one-dimensional  applications  but  cannot  be 
extended  to  two-dimensional  signals.  In  general,  a  two-dimensional  spectrum  is  not 
factorable  into  two  polynomials.  Although  progress  has  been  made  in  factoring  two- 
dimensional  spectra  in  terms  of  the  complex  cepstrum  and  the  results  are  theoretically 
exact,  these  theoretical  results  cannot  be  attained  in  practice.  We  will,  therefore,  focus 
our  attention  on  a  second  method  of  attaining  the  prediction  error  filter. 

The  solution  method  dealt  with  in  this  thesis  begins  with  the  individual  samples  of 
the  random  process  and  develops  the  prediction-error  filter  from  these  samples.  This 
method  is  known  as  autoregressive  model  fitting  or  all-pole  modeling.  The  solution  will 
be  developed  in  terms  of  orthogonal  polynomials  that  will  not  only  yield  a  Levinson-like 
algorithm  but  are  easily  extended  to  two-dimensional  signals. 


B.     ONE-DIMENSIONAL  LINEAR  PREDICTION 

1.     All-Pole  Model 

We  are  given  a  finite  length  sequence  of  the  past  values  of  the  random  sequence 
Y(n).  The  problem  is  to  estimate  Y(n)  based  on  a  linear  combination  of  these  samples. 
For  the  simplest  case  we  can  estimate  Y(n)  based  solely  on  the  previous  sample.  The 
estimate,  denoted  by  Y(n),  is 

Y{n)  =  aY(n-l)  (1) 

The  prediction  will  have  some  error,  e(n),  which  will  be  minimized  in  some  sense 
by  the  choice  of  a.  The  least-squares  minimization  criterion  is  commonly  used.  This 
criterion  provides  a  simple  approach  to  the  solution  of  the  linear-prediction  problem  and 
thus  the  modeling  problem.   We  proceed  as  follows. 

First,  we  redefine  the  estimate  as 

Y{n)  =    -aY(n-  1)    .  (2) 

Now  the  prediction  error  will  be 

e{n)  =   Y{n)  -   Y(n) 

=   Y(n)  +  aY(n-  1)    . 

The  squared  error  is  to  be  minimum  with  respect  to  the  coefficient  a.  This  is 
easily  accomplished  using  derivatives  as  shown  below.  Differentiating  [V(/?)]  with  re- 
spect to  a  gives 

—  [e2(n)l  =  2e  —e    =  0  (4) 

da  L    l  ,j  n  da    n  K  } 


or 


\Y{n)  +  aY(n  -  1)]  -f-  [7(«)  +  aY(n  -  1)]  =  0  (5) 

da 


which  vields 


lY(n)  +  aY(n-  1)][F(«-  1)]  =  0    .  (6) 

Finally,  solving  for  a 


Y{n)Y{n-\) 

a  Y{n-\)Y(n-\)      '  [) 

Since  the  value  of  a  must  minimize  the  error,  e(n),  for  all  values  of  n,  the  value 
of  a  in  (7)  could  be  rewritten  as 

n 

£V(/)r(/-i) 

-a  =    -r* (8) 

£r(/_i)r(/_i) 

Note  that  (8)  will  select  a  such  that  £e2(«)  is  minimum.   It  can  be  seen  that  the 

n 

coefficient  a  is  a  function  of  the  autocorrelation  elements  of  the  given  sample  space,  and 
this  method  is  usually  referred  to  as  the  autocorrelation  or  Yule-Walker  approach. 
If  the  order  of  the  predictor  is  now  increased  such  that 

Y(n)  =    -laJ{n-\)  +  a2Y{n-2)  +  ••■   ]  (9) 

the  prediction  error  becomes 

en  =   Y{n)-  Y(n)  =   Y{n)  +  aj{n-  1)  +  a2Y(n  -  2)  +  ■•■■    .  (10) 

If  (10)  is  considered  to  be  the  difference  equation  of  an  all-zero  filter,  the 
prediction-error  filter  can  be  described  by  taking  the  z  transform  of  (10)  as 

E(z)  =   Y(z)A{z)  (11) 

where 

A(z)  =   1  +  a,z_1  +  a2z~2  +  ....    .  (12) 

Now  if  (10)  is  rewritten  and  e(n)  is  taken  to  be  the  input  to  the  system,  the  es- 
timated signal  generator  model  results 

Y(z)  =  E(z)B(z)  =    -j^E{z)  (13) 

A(z) 

where 


B{z) 


1   +  axz       +  a2z       +  •••• 


(14) 


It  can  be  shown  that  en  is  a  white  noise  sequence  with  variance  £[<?*].  This  im- 
plies that  Y(z)  and,  consequently,  Y(n)  can  be  generated  from  a  white  noise  sequence 
as  shown  in  Fisure  1. 


Y<z) 
E(z>- 

-^  \ —  /  < 

A(zi 

*^E  (z) 

(15) 

A(z) 

=    1  +  a,z       +  a2z 

Biz) 

B(z)  = 

(16) 

1 
A{z) 

i 

!                   -1                   -2 

1  —  fl,z       —  a2z       —  ■•■ 

Figure   1.      Prediction  error  and  system  model  filters. 


This  assumes  that  the  sequence  A(z)  is  causal  and  causally  invertible.  If  this  is  the  case, 
we  can  generate  the  error  sequence  from  Y(n)  and  vice-versa  by  an  invertible  filter  op- 


eration  [Ref.  1].  The  signal  generator  model,  B(z),  is  the  best  estimator  in  a  least-squares 
sense  and  A(z)  is  the  prediction-error  filter.  These  filters  can  be  used  in  a  number  of  ways 
and  the  simple  solution  for  the  filter  coefficients  is  the  primary  reason  for  the  popularity 
of  linear  prediction  in  signal  processing. 

2.     Normal  Equations  and  Orthogonality 

Recalling  that  we  wish  to  minimize  the  sum  of  squares  of  the  error  sequence 
leads  directly  to  an  important  result  in  linear  prediction  known  as  the  orthogonality  re- 
lationship between  the  sample  space  and  the  error.  This  principle  states  that  the  error 
e(n)  is  orthogonal  to  all  the  samples  of  Y(n)  used  in  obtaining  the  estimate  Y(n).  That 
is, 

£[>(«),  Y(n  -  /)]  =  0,     1  <  i  <  N  (17) 

where  N  is  the  order  of  the  predictor.  This  result  is  derived  for  the  two-dimensional  case 
in  Chapter  II.  It  is  the  central  idea  behind  the  solution  of  the  linear  prediction  problem 
in  terms  of  orthogonal  polynomials.  It  also  leads  directly  to  the  normal  equations  in 
both  one-  and  two-dimensional  formulations. 

The  derivation  of  the  normal  equations  is  a  simple  extension  of  the 
orthogonality  relationship  (17).  Suppose  we  seek  a  kth-order  predictor  and  the  value 
of  the  prediction-error  variance.  We  will  require  (k+  1)  linearly  independent  equations 
to  solve  for  the  k  prediction  coefficients  and  the  prediction-error  variance.  These  coef- 
ficients can  be  derived  as  follows. 


-   Y{n)  -   Y[n) 

=   Y(n)  +  axY(n-l)  +  a2 Y(n  -  2)  +   +  akY(n  -  k) 


(IS) 


Substituting  (IS)  into  (17)  gives 

E\_Y{n)  +  axY{n-  1)  +  a2Y(n  -  2)  +    -  ,  }'(«-/)]  =  0      \  £  i  <,  N         (19) 

If  we  now  define  the  coefficient  a0  =  1   and  note  that  for  a  stationary  process  the 
autocorrelation,  R(x),  is  given  by  £[  ?'(/:)  7(/)]  =  Rk_,  ,  (19)  may  now  be  written  as 

v 
J^a.Yin-k^Yin-i)     =0  (20) 

k=0 


,v 

A=0 


(21) 


where  N  is  the  order  of  the  predictor  being  used  and  R  is  the  autocorrelation  function 
of  the  random  sequence  V(n).  The  N  equations  generated  by  (21)  can  be  written  in 
matrix  form  for  N  =  2  as 


R-i    RQ    /?-, 

R2      Ry      R0 


r«o 

r°i 

a, 

Ln2. 

_0_ 

(22) 


The  error  variance  c  is  given  by 


c  =  E\_e2n~]  =  E[en,  Y(n)  -   Y(n)] 
=  Elen,  Y[n)l  . 


(23) 


Substituting  (18)  into  (23)  yields 


c  =  E 


YakY(n  -  k),Y(n) 


k=0 


(24) 


which  in  terms  of  the  autocorrelation  becomes  c  =  X"*^-a  •    l:°r  example,  with  N=2, 
this  can  be  combined  with  (22)  to  yield 


Ro    R-t    R-2 
Rx     Ro    R-\ 


R,     R< 


Ro   I 


r     n 

% 

€ 

«l 

= 

0 

L«2_ 

_0_ 

(25) 


This  autocorrelation  matrix  is  seen  to  be  (N+  1)  by  (N+  1),  Toeplitz,  and  since 
R_k  =  Rk,  it  is  symmetric.  The  various  solutions  for  (25)  occupy  much  of  the  literature 
on  linear  prediction.  For  our  purposes,  the  desired  solution  must  provide  a  simple  ex- 
tension to  two-dimensional  signals  and  be  adaptable  to  a  lattice  structure.  J.H.  Justice 
[Ref.  2]  has  derived  such  a  solution  in  terms  of  Szego  polynomials.  Although  quite 
similar  to  the  Levinson  algorithm,  it  has  the  important  advantage  of  not  relying  on  the 
symmetry  of  the  correlation  matrix.  The  two-dimensional  autocorrelation  matrix  will 
be  found  to  be  Toeplitz  in  blocks,  and,  although  this  provides  some  simplification,  the 


symmetry   of  (25)  is  lost.      Fortunately,   a   solution  in  terms   of  orthogonal  Szego 
polynomials  is  easily  derived  [Ref.  3]. 

3.     Applications  of  the  Autoregressive  Model 

Linear  prediction  plays  an  integral  role  in  the  solution  of  a  number  of  engi- 
neering problems.  The  resulting  model  can  be  used  for  high  resolution  spectral  esti- 
mation [Ref.  4].  Speech  signals  can  be  reduced  from  a  large  volume  of  data  to  a 
comparatively  small  data  space  by  linear  predictive  coding  [Ref.  5].  A  model  of  the 
generating  function  for  the  random  sequence  can  be  derived.  Lattice  filters  can  be  im- 
plemented which  can  be  used  to  form  the  prediction-error  filter  in  the  forward  direction 
or  synthesize  the  random  process  in  the  inverse  direction. 

The  many  applications  of  linear  prediction  in  time  series  analysis  has  generated 
a  great  deal  of  interest  in  extending  the  results  to  multi-dimensional  random  processes. 
Multi-dimensional  signals  occur  in  situations  where  the  samples  are  spatially  dependent 
such  as  radar,  sonar,  and  image  processing.  The  signals  considered  in  this  thesis  are  wide 
sense  stationary  and  scalar-valued  sequences  of  two  variables.  The  two-dimensional 
autoregressive  model  that  is  derived  will  be  used  in  some  of  the  same  applications  clas- 
sically associated  with  the  1-D  model  including  spectral  estimation  and  sequence  gener- 
ator estimation.  An  alternative  implementation  of  the  model  in  a  lattice  structure  will 
also  be  presented. 

C.     NOTATIONAL  CONVENTIONS 

1.  Matrix  Notation 

This  thesis  will  deal  with  a  finite  size  rectangular  sample  space.  The  individual 
samples  will  be  denoted  by  Y(nx  —k,  n2  —I),  where  k  and  /  are  measures  of  distance  from 
the  point  at  which  the  value  of  the  sequence  is  being  estimated.  The  sample  space  will 
contain  NT  x  X2-1  samples.  Figure  2  shows  an  example  of  the  sample  space  with 
\1  =  \2  =  4.  In  this  case  we  would  be  predicting  the  value  of  Y(3,3)  based  on  the  15 
samples  available.  Notice  that,  in  this  case,  k  and  /  have  a  maximum  value  of  3.  Spa- 
tial samples  are  indexed  by  nl  and  n2  to  denote  rows  and  columns,  or  the  vertical  and 
horizontal  directions,  respectively. 

2.  Causal  Samples 

In  one-dimensional  signals  the  clear  meaning  of  past  and  future  leads  to  a 
straightforward  definition  of  causality.  In  two-dimensional  signals  the  ordering  of  the 
points  is  of  central  importance;  therefore  a  reference  for  past  and  future  needs  to  be  es- 
tablished.   In  this  thesis  causal  samples  will  be  taken  to  be  in  the  quarter  plane  shown 
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Figure  2.      Sample  space  for  N 1  =  N2  =  A. 


in  Figure  3.  This  is  equivalent  to  describing  the  desired  filter  as  a  third-quadrant  filter. 
Due  to  the  symmetry  of  the  2-D  autocorrelation  the  third  and  first  quadrants  produce 
the  same  results.  However,  causal  will  be  used  only  in  reference  to  the  samples  shown 
in  Figure  3. 

3.     Polynomial  Notation 

Two-dimensional  polynomials  will  be  denoted  by  Pkl  where  k  is  the  degree  of  the 
polynomial  in  the  nl  direction  and  /  is  the  degree  in  the  n2  direction.  When  the  need  to 
refer  to  an  individual  clement  in  a  polynomial  arises,  the  clement  will  be  identified  as 
Pkl'J),  where  i  and  j  give  the  power  of  the  bivariatc  variable.  The  correspondence  be- 
tween the  z  transform  of  the  system  model  and  the  coefficients  in  the  orthogonal 
polynomials  leads  to  the  use  of  z,  and  z2  as  the  independent  variables  in  this  study. 
Thus  the  polynomial  P22  will  be  a  3x3  bivariate  polynomial  consisting  of  the  elements 
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Figure  3.      Caus.il  Data  Space 


P22  = 


\  ^22(0,0) 
W2,0)zf     />22(2,1); 


/'22(°.1)-2  /,22(°»2)z2 

ft2(l,l)z,z2     /722(l,2)r,r2 

lz2      /722(2»2)71Z2 


(26) 


The  development  of  Pkl  is  in  terms  of  positive  powers  of  z.  The  desired  filter 
must  be  causal  and  therefore  in  negative  powers  of  z.  This  is  easily  accomplished  by 
multiplying  Pkl  by  z\ *zj'.  The  resulting  causal  filter  will  be  denoted  Ak,  .  The  leading 
coefficient  of  PA/  is  pk,(k,f)  which  will  become  akl(0,0)  in  the  causal  filter. 


II.     TWO-DIMENSIONAL  LINEAR  PREDICTION 

A.     PROBLEM  DEFINITION 

The  two-dimensional  linear  prediction  problem  is  conceptually  identical  to  the  one- 
dimensional  problem.  That  is,  we  attempt  to  predict  the  value  of  Y(nl,n2)  based  on  a 
linear  combination  of  the  previous  NlxN2  rectangular  block  of  data.  The  coefficients 
which  result  in  the  optimum  prediction  will  comprise  a  linear,  shift  invariant,  causal, 
prediction  error  filter.  The  problem  is  to  find  the  coefficients  that  minimize  in  some 
sense  the  error  between  the  actual  value  of  Y(nl,n2)  and  the  predicted  value,  T(nl,n2). 
This  will  require  the  solution  of  a  block-Toeplitz  matrix  which  contains  NlxNl  blocks 
and  each  block  will  contain  N2xN2  correlation  coefficients.  We  will  require  the  solution 
to  provide  the  means  for  a  direct  realization  of  the  filter  or  implementation  in  a  lattice 
structure.  The  mean  square  error  criterion  is  again  utilized  to  obtain  the  optimum  pre- 
diction filter.  If  a  solution  to  the  desired  filter  can  be  achieved  recursively,  then  the  re- 
lationship between  successively  higher  order  filters  can  be  used  to  implement  the  lattice 
filter. 

We  begin  by  seeking  a  solution  to  the  linear  prediction  problem,  and  will  assume 
that  V(nl.n2)  is  a  two-dimensional  stationary  random  sequence.  Given  the  causal  block 
of  NlxN2  data  samples  we  need  to  generate  T(nl,n2)  from  a  linear  combination  of  the 
data  that  minimizes  the  prediction  error.   The  prediction,  }7(nl,n2),  can  be  written  as 

Y(n\.n2)  =    -  al0Y(nl  -  1,«2)  -%M^2  -  1)  -  auY{n\  -  1,«2  -  1)  -  - 

-    -%i-i,^2-i  Y{n\  -  Nl  +  \,nl  -  N2  +  1) 

where  the  asl's  are  the  filter  coefficients  associated  with  data  samples  that  are  delayed  by 
5  in  the  nl  direction  and  t  in  the  n2  direction.  The  coefficients  can  be  assembled  in  a 
matrix  to  yield  the  desired  two-dimensional  linear  prediction  filter.  This  filter  will  be  an 
all-zero  model  with  the  finite  sample  space  Y  as  the  input  and  the  prediction,  Y,  as  the 
output  as  described  below. 
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Figure  4.      Generation  of  the  prediction  error 


The  prediction  process  can  be  depicted  as  in  Figure  4.    The  polynomial  A(z,  ,  z2)  is 
simply  the  z  transform  of  (27)  given  by 


A(z{,z2)  =    ~{al0z{  '   +  a0]z2  '  +  a,,z,  'z2  2  +    ••••)    . 


(29) 
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For  any  choice  of  the  coefficients  some  error  will  be  present.  This  prediction  error  is  as 
seen  in  Fisure  4  and  can  be  written  as 


AT1-1  #2-1 


e=Y{nl,n2)  +    ^    ]T  astY{nl  - s,n2 -  t)  .  (v)  ^  (0,0)  (30) 

s=0      f=0 

We  will  define  the  filter  coefficient  a00  =  1  and  condense  (30)  to  the  form 

m-i  N2-i 
e  -    £    £axrr(«l-v*2-r)    .  (31) 

5=0        f=0 

The  task  now  is  to  find  the  filter  coefficients  that  minimize  the  error. 

B.     TWO-DIMENSIONAL  ORTHOGONALITY 

1.  Description 

A 

The  error,  Y(nl,n2)-y(«l,«2),  can  be  visualized  as  in  Figure  5.  Any  linear 
combination  of  Y(nl-l,n2)  and  Y(nl,n2-1)  will  yield  an  estimate  of  Y  which  lies  in  the 
plane  of  Y(nl-l,n2)  and  Y(nl,n2-1).  If  Y(nl,n2)  lies  in  this  plane,  then  Y  will  be  exact 
and  there  will  be  no  error.  For  any  other  Y(nl,n2)  the  error  will  be  minimum  if  it  is 
orthogonal  to  the  sample  space  spanned  by  Y(nl-l,n2)  and  Y(nl,n2-1).  This  is  equiv- 
alent to  defining  Y  as  that  portion  of  Y(nl,n2)  that  is  parallel  to  the  sample  space.  Thus 
Y-Y  leaves  only  the  component  of  Y  that  is  perpendicular  to  the  sample  space  as  the 
error.  Although  it  can  not  be  drawn  for  higher  order  predictors,  the  minimum  error  is 
always  such  that  it  is  orthogonal  to  the  space  spanned  by  all  the  samples  of  Y  used  in 
the  estimate. 

2.  Derivation 

The  orthogonality  of  the  error  to  the  sample  space  can  be  easily  derived.  First 
we  define  the  inner  product,  <X,Y>  ,  of  two  sequences,  X  and  Y,  by 

<x,y>=y,Y,x^-  w 

k       I 

Now  the  inner  product  of  the  squared  error  can  be  expressed  as 

<e2>  =  <(Y-  y)2>  (33) 

and  we  seek  a  minimum  squared  error  with  respect  to  the  filter  coefficients  a. 
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Figure  5.      A  visualization  of  the  orthogonality  principle  for  two  samples. 


-4r<c2>=0 
da 


(34) 


or 


<2e-^-e>  =  0 
da 


(35) 


where  — —  are  the  partial  derivatives  with  respect  to  the  elements  in  the  coefficient  ma- 
da 

trix.    Substituting  (30)  for  e 


d 


M-l  A7-1 


e  =  -i-lY{n\tn2)+  Y,    X*„n«l-*,«2-/)] 


da  da 


(36) 


5=0     r=0 


leads  to 
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-4r  e  =  [Y{n\  -  \,n2)J{n\,n2  -  1),  ••■•  ,  Y{n\  -  AT  -  \,nl  -  N2  -  1)]  .         (37) 
da 

Substituting  (37)  into  (35)  yields 

<e,Y(n\,-s,n2-t)>    =  0  ,    0  <  s  <  AT 

0<r<A'2  (38) 

V^O  . 

This  result,  (38),  states  that  the  inner  product  of  the  error  and  each  of  the  past 
samples  of  Y(nl,n2)  used  in  the  estimate  is  zero  which  requires  that  the  error  be 
orthogonal  to  the  sample  space.  Not  only  is  the  error  orthogonal  to  the  sample  space, 
but  errors  created  by  successively  higher-order  predictors  are  orthogonal  to  each  other. 
Also,  the  sequence  e(nu  n2)  will  be  uncorrelated  to  the  degree  that  the  order  of  the  pre- 
dictor allows.  That  is,  if  we  had  access  to  an  infinite  causal  data  space,  then  we  would 
produce  a  purely  white  error  sequence.  Since  we  will  always  be  limited  to  some  finite 
size  filter,  we  will  produce  an  error  sequence  that  is  partially  correlated. 
3.     T>vo-Dimensional  Normal  Equations 

The  sequences  considered  in  this  thesis  are  taken  to  be  wide-sense  stationary. 
If  they  are  also  taken  to  be  zero-mean  signals,  the  inner  product  of  two  sequences  is 
equivalent  to  the  expected  value  of  the  product  of  the  sequences.  The  expected  value 
of  Y(nl,n2)  and  Y(nl-s,n2-t)  reduces  to  a  two-dimensional  autocorrelation  that  is  solely 
a  function  of  the  delays  s  and  t.   Thus. 

<  Y{nl,n2),Y{n\  -  s,n2  -  t)>    =  E[Y{n\,n2),Y{n\  -  s,n2  -  /)]  =  Rs!  (39) 

where  R  will  be  used  to  denote  the  autocorrelation  with  lags  s  in  the  nl  direction  and  t 
in  the  n2  direction.  Equation  (38)  can  now  be  used  to  form  a  block-Toeplitz  matrix  that 
is  analogous  to  the  Toeplitz  matrix  generated  by  the  normal  equations  in  one- 
dimensional  linear  prediction.    Thus,  we  can  write 

E[e,  Y{n\  -  s,n2  -  /)]  =  0     (s.t)  *  (0,0)  (40) 


substituting  for  e 

E 

/Vl-l  N2-1 

Y(nl,n2)+  ]T    £ 
fc=0     /=o 

aklY(n\  -  k,n2  -  I) ,  Y{n\  -  s,n2  -  /) 


=  0    0<s<AT 

0<t<N2     (41) 
(s,t)  #  (0,0) 
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that  can  be  written  as 


A'l-1  A"2-1 


0  <  5  <  /Vl 

y    ^  aklR(s  -k,t-l)  =  0    .  0  <  /  <  A'2 


(42) 


A-=0      /=0 


For  example,  with  a  N1  =  N2=  2,  and  with  the  first  filter  coefficient  defined  as 
a00  =  1  the  normal  equations  become 


E\  Y{n\,n2)  +  a0l  Y(nl,n2  -  1)  +  a]0Y(n\  -  1,«2)  + 


+  auY{n\  -  1,«2  -  1) ,  Y{n\  -  s,n2  -  t)\   =  0 


(43) 


where  (42)  or  (43)  are  the  two  dimensional  normal  equations.  They  are  a  direct  result 
of  the  orthogonality  relationship  described  by  (40).  These  equations  can  be  written  out 
to  yield  NlxN2-l  linearly  independent  equations.  For  example,  with  N1  =  N2=2,  we 
have 


s  =  0j=  1      /?01  +a0lR00  +  a]0R_]A  +  anR_h0  =  0 


(44) 


5  =  1,/  =  0       Rw  +  fl0l^l,_l  +  «io^oo  +  auRo-\    =   ° 


(45) 


s=\j=\         Ru  +a0]R]0  +  awR0]  +  anRQQ  =  0 


(46) 


Equations  (44)-(46)  can  be  augmented  by  the  expression  for  the  squared  error 
c  which  is  eiven  bv 


e  =  £|V]  =  £[e,  Y-  Y] 


(47) 


or 


=  E 


A']-l  A'2- 


k=0     1=0 


V  aklY(n\  -  k,n2  -  !),Y{n\,n2) 


(48) 


and  in  terms  of  the  autocorellation 


A'l-l  A/2-1 
k=0     1=0 


Assembling  the  above  equations  into  a  matrix  yields 


(49) 
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#00     #0,-1 

#-1,0 

#-l,-l~ 

a00 

€ 

#01        #00 

#-1,1 

#-1,0 

«01 

0 

#10     #1,-1 

#00 

#0,-1 

a\Q 

0 

#11        #10 

#01 

#00 

_an_ 

_0_ 

(50) 


The  resulting  matrix  is  (NlN2)x(NlN2)  square  and  will  consist  of  blocks  which 
are  Toeplitz.  There  will  be  NlxNT  blocks  and  each  block  will  contain  N2xN2  corre- 
lation coefficients.  The  matrix  can  be  seen  to  be  Toeplitz  with  respect  to  the  blocks  and 
can  be  written  as 


c£0     O 


//■ 


o 


% 


(51) 


where  the  <£>k's  are  blocks  in  the  autocorrelation  matrix  and  H  indicates  hermitian 
transposition.  Each  O*  is  an  N2xN2  Toeplitz  matrix.  For  the  example  from  (50),  we 
have 


(52) 


(53) 

(54) 
(55) 

(56) 

If  an  element  in  <t>k  is  Rst  ,  the  corresponding  element  in  <&HisR'_sr.  Since  we 
are  dealing  only  with  real  data  the  conjugate  notation  can  be  dropped  from  R.  It  is  also 
worthwhile  to  note  that   Rsr  =  R_s_,  . 

The  solution  of  this  matrix  for  the  coefficients  a!t  will  yield  the  two-dimensional 
whitening  filter.  There  has  been  a  great  deal  of  research  involving  block  Toeplitz  ma- 
trices and  there  are  several  algorithms  that  could  be  used  to  obtain  the  coefficients.  In 
order  to  implement  the  filter  in  a  lattice  structure,  a  recursive  solution  beginning  with  a 
zeroth-order  filter  and  proceeding  to  the  desired  order  is  necessary.    In  one-dimensional 


D    — 

"#00 
_#01 

#0,-1 
#00 

1    = 

#10 

_#11 

#.,-1 

#10 

a0    =    («00, 

a\  =  (a\o, 

«ll) 

c   ■■ 

=  (£, 

o)r 

16 


prediction  error  filtering,  the  Levinson  algorithm  successively  represents  the  higher  order 
filters  in  terms  of  reflection  coefficients.  It  will  be  seen  that  an  analogous  solution  to  the 
two-dimensional  prediction  error  filter  can  be  developed  in  terms  of  Szego  polynomials. 

C.     SZEGO  POLYNOMIALS 

The  solution  to  the  system  of  equations  Oa  =  c  does  not  require  that  the  inverse  of 
the  matrix  O  be  found.  All  that  is  required  is  an  orthogonal  set  of  polynomials  that  span 
the  space  defined  by  O.  This  result  follows  the  orthogonality  relationship  of  the  pre- 
diction error  and  the  sample  space.  For  example,  if  we  had  a  second  order  and  a  first 
order  prediction  error  filter  the  corresponding  prediction  errors  would  be 

ex  =   Y(\)-ax(\)Y(0)  (57) 

e2  =   Y(2)-a2(l)Y(\)-a2(2)Y(0)  (58) 

EL>2.«i]  =  Efo.ni)-  «i(OT0)]  (59) 

But  (44)  requires  that  e2  is  orthogonal  to  both  Y(l)  and  Y(0),  and  therefore 

£LVi3  =  o  .  (60) 

It  is  evident  from  the  above  observations  that  the  choice  of  AA_,  must  be  such  that 
it  is  orthogonal  to  AA.  If  we  begin  by  defining  A0  as  any  convenient  value  and  proceed 
to  find  the  successively  higher-order  orthogonal  polynomials  until  the  desired  predictor 
order  is  reached,  the  final  polynomial  will  hold  the  coefficients  which  will  provide  the 
optimum  prediction-error  filter.  To  accomplish  this  we  will  use  Szego  polynomials. 
Szego  polynomials  are  a  special  class  of  polynomials  which  have  zero  inner  products 
over  an  interval  on  the  unit  circle  [Ref.  6].  This  is  precisely  the  requirement  of  the 
orthogonal  polynomials  relating  successively  higher  order-prediction  errors.  J.H.  Justice 
[Ref.  3]  has  presented  a  recursive  method  of  generating  Szego  polynomials  of  the  form 

pk  =  P{k)zk  +  Pk-\z>C~l  +  -■    +Po  (6I) 

where 

z  =  dkB  (62) 

The  matrix  O  is  an  inner  product  matrix  which  may  be  utilized  to  orthonormalize  a 
sequence  of  polynomials  on  the  unit  circle  [Ref.  6].   The  resulting  orthonormal  sequence 
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will  be  composed  of  Szego  polynomials.  If  the  filter  sequence  Xk  is  determined  in  terms 
of  Szego  polynomials  that  are  orthogonal  over  the  space  associated  with  the 
autocorrelation  matrix  ,  the  prediction  error  filter  will  result  [Ref.  2].  Further,  the  Szego 
polynomials  can  be  determined  in  a  recursive  manner  which  is  necessary  for  the  repre- 
sentation of  the  filter  in  a  lattice  structure.  In  a  one-dimensional  filter  this  will  result  in 
the  Levinson  algorithm.  J.H.  Justice  [Ref.  2]  has  developed  the  bivariate  extension  of 
Szego  polynomials  that  can  be  utilized  to  orthonormalize  the  block  Toeplitz  matrix 
generated  by  the  two-dimensional  autocorrelation  of  the  sequence  Y(nl,n2). 
1.     Levinson  Algorithm  in  Terms  of  Szego  Polynomials 

Any  positive  definite  matrix  can  be  defined  as  an  inner  product  space  and  can 
thus  be  used  to  define  orthogonal  polynomials  spanning  that  space.  The  autocorrelation 
matrix  produces  such  a  space.  If  the  autocorrelation  matrix  is  generated  from  a  se- 
quence, Y=  Y0,  Yi,-,Yjf,  then  the  sequence  of  polynomials,  l,z,  z2,  ■■  ,zN  can  be 
orthonormalized  with  respect  to  the  autocorrelation  matrix.  The  resulting  polynomials 
are  P0{z),  Px{z),  •-•,  Ph{z).  The  subscript  on  P  denotes  the  degree  of  the  polynomial.  The 
requirement  on  PA  is  that  it  is  orthogonal  to  all  polynomials  of  lesser  degree.  Thus,  we 
have 

<Pk,P1>   =  S(k-[)  (63) 

A  one-dimensional  sequence  of  orthogonal  polynomials  will  consist  of  the  terms 

A)  =  Poo 

P\  =P\\2    +    P\0 

(64) 


A' 


-V-l 


Pn  =  Psnz     +  Pnn-\z  '     +•-  +Pxo 
and  in  general, 


=  YpJ    ■  (65) 


Pk  :  :  /JPkj 
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Justice's  method  for  determining  these  polynomials  recursively  results  in  an  al- 
gorithm equivalent  to  the  Levinson  algorithm.  This  algorithm  for  the  1-D  case  is  briefly 
summarized  as  follows.   We  begin  by  setting 

P0  =  #0-T    .  (66) 

Then  the  subsequent  polynomials  can  be  determined  by  the  recursive  relationship 

Pn+l{z)  =  (1  -  I  /  J  2pl)-T(zPn(z)  -  >.,fnnPrn{z))  (67) 

where 

Prn(z)  =  znPn(z-])  (68) 

and 

n 

)-n    =    YjRk+lPnk     ■  (69) 

fc=0 

The  equation  (67)  can  be  expressed  in  the  form  of  a  causal  prediction-error  filter 
that  yields  the  familiar  form  of  the  Levinson  algorithm.  The  nth-order  polynomial 
produced  by  (67)  will  be  of  the  form: 

Pn(z)    =   Pnnz"  +  Pnn-l2"'1  +  -■    +  PnO     •  (70) 

The  causal  prediction-error  filter  is  of  the  form: 

An{z)  =  an0  +  an]z~]  +■-  +  annz~n  (71) 

where  we  let  anp  -  pn/,_p  and  A„  =  i~nPn  .   Then  multiplying  (67)  by  z-("+1)  yields 

z-(n+l)Pn+](z)  =  z-"Pn(z)  -  z-<n+%Prn{Z)  (72) 

where  the  normalizing  factor,  (1  —  U„l2/>™)_7,  has  been  absorbed  into  the  definition 
of  ).n  and  the  leading  coefficients  are  defined  as  pm  =  a„0  =  1.  Now,  we  define  the 
backwards  polynomials  as  follows 

Prn(z)  =  znPn(z-])  (73) 
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Arn(z)  =  z  "An(z     ) 


=  z-\znPn(z-])) 


(74) 
Prn(z)  =  znArn(z)  (75) 


where  (72)  can  now  be  written  as 


An+](z)  =  An{z)-z-x).nArn{z)  (76) 


4+1(z)  =  z-^\An{Tl)  -  zXnArn{z~1)) 

=  2~^Arn{z)-z~n).n{znAn{z)) 


(77) 
4(z)  =  z'lArn(z)-AnAn(z)  (78) 


K  =    ^^^ (79) 

/=0 

Equations  (76),  (78)  and  (79)  represent  the  usual  Levinson  recursion  of  one- 
dimensional  linear  prediction  problems.  The  coefficients  found  using  the  Levinson  al- 
gorithm will  be  developed  in  negative  powers  of  z.  The  coefficients  found  using  (67)  are 
in  positive  z  and  differ  by  a  normalizing  factor.  Both  sets  of  polynomials,  A  and  P  ,  will 
be  orthogonal  sets  of  polynomials  on  the  unit  circle  spanning  the  same  inner  product 
space.  Since  the  definition  of  the  ?.„  is  simpler  using  (67),  we  will  proceed  to  the  two 
dimensional  extension  using  positive  powers  of  z  and  make  the  necessary  conversion  to 
a  causal  filter  after  the  coefficients  have  been  found. 
2.      2-D  Extension 

Now  we  want  to  extend  these  ideas  to  two-dimensional  filters.  The  necessary 
extension  to  bivariate  orthogonal  polynomials  has  been  developed  by  Justice  [Ref.  2]. 
The  derived  polynomials  will  be  orthogonal  in  two  variables  with  respect  to  the  inner 
product  space  of  O  as  described  by 

£[^*/.^r]  =  S(k-s)S(l-t)    .  (80) 
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The  set  of  polynomials  to  be  orthonormalized  are  shown  in  Figure  6.  The  procedure 
starts  with  the  first  row  which  is  a  function  of  z2  alone.  This  row,  therefore,  can  be 
orthonormalized  by  the  same  method  presented  for  one-dimensional  systems. 
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Figure  6. 


Variables  used  to  form  the  bivariate  orthogonal  polynomials 


This  is  equivalent  to  taking  the  block  O0  from  the  two-dimensional  autocorrelation  ma- 
trix and  finding  the  filter  coefficients  a0lt  aQ2,  •■■  ,  c^m.   Thus, 


<I>0  = 


Koo     R 


Re 


0,-1 
*00 


R 


m0 


R 


0,-m 


R 


0,-<m-l) 


R 


00 


(81) 


where 


ao  "~  (a00>  a01'  ""  '  GQm) 


(82) 


It  is  easily  shown  that  R0m  =  R0_m,  so  that  this  is  a  symmetric  matrix  and  a  re- 
cursive solution  has  been  presented.  Once  the  first  row  of  orthogonal  polynomials  has 
been  determined,  the  polynomials  of  the  next  row  can  be  expressed  as  a  linear  combi- 
nation of  the  lower-order  polynomials.  For  example,  with  N2  =  3  the  first  row  will  have 
polynomials  PM  ,  P0l  ,  and  P02.  The  first  polynomial  in  the  second  row  to  be  found  is 
/>,„.   This  can  be  accomplished  as  follows: 
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First  note  that  if  P00  is  multiplied  by  z,  its  order  is  changed  and  it  is  no  longer  a 
member  of  the  orthogonal  set.   We  may  then  write 

Z1M)0    =    ^-00  Ao    +    *01°01    +    ^02-^02  (83) 

and  after  multiplying  on  both  sides  by  P^  and  taking  the  expected  value  we  find 

£[z,?oo.  Pld    =    ^OO^C^O^IO]    +    ^0l£C^01.^10]    +    ^02^[^02^]0]     •  (84) 

Since  P]0  is  orthogonal  to  all  polynomials  of  lower  order  in  zu  this  reduces  to 

E^z^Pqq.  P10]  =  A00  (85) 


and  similarly 


fE^oo^oi]  =  ;-oi  (86) 

£[^1^00.^02]  =  ^02    •  (87) 


Thus,  in  general  we  find  that 


/-l  k     A7-1 


z\pki  =  }~klPk+i,l  +  /^A'^+U  +  lu  )_^^-stPsf  (88) 


:=0  s=0    1=0 


Now  we  define  the  reverse  polynomial 

0<j<m  (89) 

m  =  N2  -  1 

Note  that  P\,  will  also  be  an  orthogonal  set  of  polynomials  and  will  span  the  same  space 
as  P,j.   Now  (SS)  may  be  written  in  terms  of  the  reverse  polynomials 

/- 1  A-      m 


t=0  5=0  /=0 


where 


But  for    0  <  5  <  k  -  1 


A'st=  ElZ]PkhPrst-]    .  (91) 
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£[>  A, />;,]=  0  (92) 

so  that  (90)  may  be  reduced  to 


i-\ 


z\Pki~  ^ki^k+\,i  +  /^^kt^k+\,t  +  /  A'kt^kt  (93) 


t=0  t=0 


where 


and 


4,=  EizxPkbPk^  (94) 

>'k<=  ELziPuK']    ■  (95) 

The  polynomials,  Pk+U ,  may  now  be  found  from  the  lower  order  polynomials 


/ —  1  rn 

Pk+\,l  =   z\Pkl  ~    /  /■kt^k+\,t  ~   /_J-'ktPkt  (96) 


;=0  i=0 


W/-kfk+\,l\\ 


This  completes  the  development  of  the  bivariate  Szego  polynomials  as  pre- 
sented by  Justice  [Ref.  2].  The  polynomial  PH  will  hold  the  coefficients  for  the  k  xl 
prediction  error  filter.  The  desired  causal  prediction  error  filter,  Aw,  is  obtained  by 
multiplying  Pk!  by  z\kz^  and  dividing  by  pk,{k,[)  to  force  the  leading  coefficient 
fl«/(0,0)  =  1  as  was  assumed.  The  signal  generator  model  is  found  from  (30)  with  e(wj«2) 
taken  as  the  input. 

B{2X,Z2)=—^   = 1 — (98) 

^K^,z2)  l+a10Zi     +a0]z2    +auz]    z2    +  .... 

This  again  assumes  that  A(zuz2)  is  a  causally  invertible  sequence. 
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D.     LATTICE  IMPLEMENTATION 

1.     Analysis  Lattice 

It  is  now  a  simple  matter  to  utilize  the  orthogonal  polynomials  for  lattice  im- 
plementation. The  first  row  is  a  function  of  z2.  only  and  the  lattice  structure  is  identical 
to  the  one-dimensional  case.  All  subsequent  polynomials  have  a  connection  with  each 
of  the  polynomials  in  the  previous  row  and  each  of  the  polynomials  in  the  same  row  but 
of  lower  order  in  z2. 

Using  the  familiar  form  of  the  Levinson  algorithm  the  first  row  of  the  lattice 
filter  may  be  drawn  as  shown  in  Figure  7,  given  by 

^0,/c+l    =    A0k   ~   z2    ^k^Ok  (99) 

and 

A0J<+\    =    z2_I4a-    -    >-kA0k     ■  (100) 

The  orthogonal  polynomials  that  have  been  developed  will  hold  the  filter  coefficients  in 
descending  powers  of  2,  and  z2.  That  is,  the  lowest  order  term  in  Pk,  is  the  term  associated 
with  the  greatest  delay  in  the  filter  difference  equation.  We  could  use  (94),  (95)  and  (96) 
to  implement  a  first-quadrant  filter  that  would  be  statistically  equivalent  to  the  third- 
quadrant  causal  filter  we  desire.  If  all  the  data  are  stored  ahead  of  time,  which  is  gen- 
erally the  case  for  2-D  sequences,  this  would  not  cause  any  computational  difficulties. 
We  will  present  the  causal  implementation  of  the  lattice  structure.  The  first  quadrant 
development  is  similar. 

In  order  to  correlate  the  powers  of  Z]andz2  with  the  delay  and  thus  have  a 
causal  filter,  we  will  define  Ak,  =  z\kz2'Pkl.  This  will  result  in  a  realizable  design  and  the 
form  of  the  lattice  filter.  The  first  row  of  the  filter  is  given  by  (99)  and  (100).  We  now 
assume  that  row  k  has  been  constructed  and  that  we  proceed  to  row  k+1.  Multiplying 
(96)  by  zf'^1^-' yields 

l-l  m 

Akl2\  z2   M:+l,/  —   zl    z2   rkl         ^//-/^l  z2   J  k+\,:         zl  z2    /  /  kt1  kt  liU1J 

that  can  be  written  in  the  form 
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Figure  7. 


First  row  of  the  analysis  lattice 
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'■kl^k+ij  —  Ak, 


l-\ 


t=o 


-('-') 


^kAk+\,t 


z\]/_j}/K'Akt 


t=0 


(102) 


Now  with  ArMJ   defined  as  Zi{k+i)z;mAk+1J(zy\  i\ '),  we  have 


i-\ 


^ki^k+u 


-1  ^r 


A-/ 


z2      /*r/^/c+i,r 


/=0 


E 

f=0 


^,/l 


/cr^/cr 


(103) 


Equations  (102)  and  (103)  provide  the  means  for  implementing  the  two-dimensional 
lattice  structure.  Before  proceeding  to  a  schematic  representation  of  the  structure,  it  is 
worthwhile  to  explore  the  relationship  of  the  Akl's  and  the  filter  coefficients.  This  will 
provide  some  insight  into  the  operation  of  the  lattice  filter  and  make  it  apparent  that  the 
/'s  are  sufficient  to  represent  the  filter.  Assume  that  we  have  already  found  the  lower 
order  filters,  Am%  Am,  A02,  and  Aw  and  now  seek  /!,,.  The  coefficients  that  are  available 
are  shown  in  Figure  8.   By  (102)  we  find  that 
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Figure  8.      The  lower  order  polynomials  for  /4, 
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^•oi^n  —  ^01  —  z2  ^-oo^io  —  zi   ^-'oo^oo- 2i   ^'oi^oi- zi   ^-'02^02    •  (104) 

The  filter  Au(zu  z2)  will  then  be  of  the  form 

An(zltz2)  =  a00  +  amz~l  +  a]0Z2l  +  anz~]  Z21  (105) 

and  (104)  can  now  be  used  to  equate  the  coefficients  in  An  with  the  lower  order 
polynomials.  This  is  shown  in  Figure  9.  The  equations  for  the  coefficients  can  be  writ- 
ten directly  from  Figure  9  as 

a,  ,(00)  =  M00)  (106) 

01,(01)  =  00,(01)  -  ;.0oMOO)  (107) 

*n(10)  =    -  >-'o2«02(02)  (108) 

«,,(!!)  =    -C/oo«,o(10)  +  >-'oi«oi(01)  +  /'o2«o2(01)]  (109) 

0  =    -  [;.'oo«oo(00)  +  /'oi«oi(00)  +  ;.'o2«o2(00)]    .  (HO) 

It  can  be  seen  that  there  are  enough  equations  to  determine  all  the  filter  coefficients  of 
An.  Conversely,  if  the  filters  are  known,  then  the  A's  can  be  completely  specified.  This 
relationship  allows  the  filter  to  be  realized  in  either  the  direct  form  using  the  filter  coef- 
ficients or  in  the  lattice  structure  using  the  X's.  If  the  same  procedure  is  carried  out 
starting  with  (102).  the  coefficients  in  A'n  are  found  to  correspond  to  those  already  found 
for  Au.   This  is  of  course  a  requirement  if  the  algorithm  is  correct. 

The  second  row  of  the  lattice  structure  described  by  (102)  and  (103)  is  shown  in 
Figure  10  for  a  2x3  filter.  The  first  row  of  the  filter  is  shown  in  Figure  7.  All  subsequent 
rows  of  the  filter  will  have  the  same  configuration.  The  structure  can  be  seen  to  be  quite 
similar  to  a  one-dimensional  lattice.  There  are.  however,  two  important  differences. 
First,  there  is  a  connection  to  all  the  lower-order  predictors  in  the  row.  Also  the  inverse 
delay  in  the  reverse  error  path  appears  to  be  a  non-causal  operator.  These  inverse  delays 
result  from  the  definition  of  Arkl=  z^z^A^z^,  z^1)  and  force  a  specific  ordering  of  points 
being  inputed  to  the  lattice.  For  example,  if  we  were  computing  the  reverse  error, 
eru(nl,n2),  for  a  filter  with  three  horizontal  filter  coefficients,  m=  2,  then  the  z  transforms 
of  of  these  errors  are  given  by 

£[,(rl,;2)  =   Y{z„z2)ATn{z,z2)  (111) 
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£[,(zl,z2)  =   Y(z],z2)(a00z]  !z2  2  +  aQ]z^  Xz2  1  +  awz2  2  +  anz2  ')  (112) 

£[0(zl,z2)  =   r(zi,z2)(a00z^1z2"2  +  a10z2"2)    .  (113) 

By  referring  to  Figure  9  we  see  that  the  coefficients  an(0,l)  and  alo(0,0)  must  be 
in  the  same  power  of  z,z2  and  likewise  au(l,l)  and  a10(l,0)  .  Therefore,  Aru{zu  z2)  must  be 
multiplied  by  z2  This  would  cause  no  problems  if  the  direct  form  implementation  was 
being  used  since 

z2£[0(zl,z2)  =   r(z1,z2)(a00z1~1z2"1 +a10z2_1)    .  (114) 

There  are  no  non-causal  terms.  In  the  lattice  structure  instead  of  linking  actual  filters 
we  link  error  outputs.  The  calculation  of  ern{n\,n2)  requires  e[0(n\,n2  +  1)  .  This  requires 
the  next  horizontal  value  of  el0  before  the  present  value  of  en  can  be  computed.  The  end 
result  of  this  is  that  the  previous  row  errors  must  be  completely  calculated  and  stored 
before  the  present  row  errors  can  be  computed.  This  forces  the  input  to  proceed  across 
rows. 

The  filter  shown  in  Figure  10  represents  extending  the  filter  order  by  rows.  If 
at  some  point  we  wish  to  extend  the  filter  by  columns  the  procedure  is  identical  and  the 
next  column  would  have  the  same  connections  as  the  row  extension  with  the  subscripts 
transposed. 

2.     Synthesis  Lattice 

The  all-pole  model  provides  the  means  to  synthesize  the  original  signal  from  the 
error  sequence.  The  synthesis  filter  can  also  be  implemented  in  a  lattice  structure. 
Multiplying  the  recursion  relationships  (102)  and  (103)  by  Y(z)  we  arrive  at  the  recursion 
for  the  error  sequences. 


^il~%tY(zl,z2)Ak+u  -  z;xYJ)'krY^z2)Ark!  (115) 

1=0  t= 


ek+\,I  =  ekl  - 


Z2{l~t)htek+\,t  -  zT1Y.;-'kAt  (116) 
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Figure  9.      Coefficient  relationship  between  orthogonal  polynomials. 


/-i 


m 


ek+\j  "  zi  **/ 


z2       ^A+I.f 


-  2_/'*'e*' 


/=0 


/=0 


(117) 


Since  Y(nl,n2)  corresponds  to  e00,  we  must  write  (116)  in  decreasing  order  re- 
sulting in 


z2(l"l)^Jaek+\,t   +   ^Zj-'kfikt 

,=0 
1=0 


(118) 


The  first  row  of  this  filter  is  again  a  function  of  one  variable  and  is  given  by 

ekl  =    ek,l+\    +    *leklz2 


(119) 


or 


ek,l+ 1  -  eki  ~  ^ieki 


(120) 
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Figure   10.      One  cell  of  a  2-D  analysis  lattice 
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This  is  realized  as  shown  in  Figure  11,  and  the  subsequent  rows  are  given  by 
(118).  The  signal  flow  graph  is  similar  to  Figure  10  with  the  signal  flow  inverted  and  the 
sign  changes  as  indicated  by  (118).  The  node  at  e00  is  shown  in  Figure  12.  Computer 
simulation  of  these  filters  and  the  results  obtained  are  presented  in  Chapter  III. 
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Figure   11.        First  row  of  the  synthesis  lattice 
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33 


III.     COMPUTER  IMPLEMENTATION  AND  RESULTS 

A.     IMPLEMENTATION 

1.  Autocorrelation  Function 

The  autocorrelation  function  referred  to  throughout  this  thesis  is  actually  an 
estimate  based  on  a  finite  sample  space.  There  are  several  ways  to  compute  the  esti- 
mated autocorrelation  coefficients.  For  a  large,  complex-valued  sample  space  the  most 
efficient  method  is  the  use  of  Fourier  transforms.  For  the  small  real-valued  data  space 
used  in  the  examples  of  this  study,  the  autocorrelation  is  estimated  directly  from  the 
available  data.   The  autocorrelation  coefficients  are  given  by 

R{n\,n2)  =  Y^jY{k,[)Y{n\+k,n2  +  [)    .  (121) 

k       I 

2.  Reflection  Coefficients 

The  basic  algorithms  for  calculating  the  coefficients  were  presented  by  Justice 
[Ref.  2].  With  some  modification  these  algorithms  are  presented  in  this  and  the  follow- 
ing section. 

The  coefficients  that  link  succesive  order  filters  are  derived  in  Chapter  II.  They 
will  be  referred  to  as  reflection  coefficients  even  though  this  term  is  descriptive  of  only 
the  first  row.   The  first  row  reflection  coefficients  were  found  to  be 

h~  /  JR0,-(k+\)P0I,0k     •  O22) 

k=0 

This  result  is  simply  the  one  dimensional  Levinson  algorithm  appended  with  a  second 
subscript. 

For  subsequent  rows  there  are  two  reflection  coefficients,  one  for  the  forward 
polynomials  and  one  for  the  reverse.   They  are  obtained  as 

Ak!=  <zlPkhPk+u>  (123) 

and 

l'ta-<*iPu.Plj>    ■  (124) 
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These  are  inner  products  over  the  autocorrelation  space  and  are  computed  as 

/  *+]     r 

hi  =  Y/kikJ%Yj  YjR^-rJ-sPk+ht(r^  (125) 


and 


m  k      I 

j=0  r=0  5=0 

where  any  negative  value  in  a  subscript  results  in  a  zero. 
3.     Filter  Coefficients 

The  first  row  of  prediction  error  filters  are  polynomials  of  the  form  P0k  .    They 
can  be  calculated  by  the  relationships 

/>o,o(0,0)  =  (y-T  (127) 

/WO,/)- (I  "  U/lVo^Or^W0*/'-  V-'-tPoJPMJLW-J))        (12S) 

0 </</+ 1 

(129) 
0  <  /  <  m  -  1 

where  m  is  the  maximum  filter  order  in  the  n2  direction  and  /  is  the  filter  order  presently 
being  computed. 

All  subsequent  rows  will  increase  the  filter  order  in  the  nl  direction  and  can  be 
recursively  computed  by 

l-\  m 

pk+ifij) = pd'  -  iJ) -  2_JiktPk+i,t(*j)  -  2^  >'kiPkAk  -  i>m  -J)  •      (13°) 

t=j  i=m-j 

0  <  /  <  k  +  1 

(131) 
0</</ 

This  recursion  can  be  carried  on  until  the  desired  filter  order  is  reached.     The  only 
polynomials  required  to  increase  the  filter  size  are  the  preceding  row  and  the  lower  order 


35 


filters  in  the  same  row.    It  may  be  desired  to  normalize  the  filter  equation,  and  this  is 
easily  carried  out  by 

Pk+\pJ)  =  —j (132) 


'k+\,l 


where 


<W=I|/W  (133) 


or 


^ir^iHli^^,.^,^)]!   .  (134) 

\/=o  j=0  / 

B.     RESULTS 

1.  System  Model 

The  algorithm  which  has  been  developed  in  this  study  is  particularly  useful  in 
finding  the  autoregressive  or  all-pole  model  of  the  system  which  has  generated  these 
data.  Several  examples  are  presented  in  this  section.  Notice  that  the  computer  gener- 
ated results  have  the  negative  of  the  desired  values.  This  results  from  our  original  defi- 
nition of  the  prediction  of  Y  in  (9)  as 

Y=-ta0]Y(n\,n2-l)  +  alQY(nl-l,n2)---l    .  (135) 

Figures  13  and  14  show  various  size  system  models  that  were  derived  from  the 
data  generated  by  the  difference  equation 

Y(nl,n2)  =  AY{nl  -  1,«2)  +  .6Y{nl,n2  -  1)  -  .6Y(n\  -  \.nl  -  1)  +  X{n\,nl) .  (136) 

Figure  13  contains  the  system  models  resulting  from  driving  the  difference  equation  with 
an  unit  impulse.  Figure  14  is  the  result  of  driving  the  system  with  a  pseudo  white  noise 
sequence.  The  discrepancies  in  the  filter  coefficients  can  be  mainly  attributed  to  insuf- 
ficient 'whiteness'  in  the  noise  data. 

2.  Spectral  Estimation 

If  we  assume  the  signal,  Y(nl,n2),  is  generated  by  the  model  of  Figure  15,  the 
power  spectrum  is  given  by 
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FILTER  ORDER  —  N1  =  2,N2  =  1 
1.0000       -.6000 
-.4000       0.6000 
0.0000       0.0000 


FILTER  ORDER  — N1  =  2,N2  =  2 

1.0000       -.6000  -.0001 

-.4000       0.6000  -.0000 

0.0000       0.0000  0.0001 


FILTER  ORDER  — N1  =  2,N2  =  3 

1-JS22       --6000       0.0000  0.0000 

"•4000       0.6000       0.0000  0  0000 

o.oooo     o.oooo     o.oooi  Sloooo 


FILTER  ORDER  — N1  =  3,N2  =  0 

1  .0000 

-.0917 

0.0000 

0.0000 


FILTER  ORDER  — N1  =  3,N2  =  1 
1.0000       -.6000 
-.4000       0.6000 
-.0000       0.0000 
0.0000       -.0000 


FILTER  ORDER  — N1  =  3,N2  =  2 
1.0000       -.6000       -.0001 
-•4000       0.6000       -.0000 
-•0000       0.0000       0.0001 
0-0000       -.0000       -.0000 


FILTER  ORDER  — N1  =  3,N2  =  3 
1.0000       -.6000 
-.4000       0.6000 
0.0000       0.0000 
0.0000      -.0000 


0 

0000 

0 

0000 

0 

0000 

0 

0000 

0 

0001 

0 

0000 

0 

0000 

0 

0000 

Figure  13.      System  models  for  Eq.136  with  an  impulse  input 
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FILTER 

ORDER  —  N1  =  2,N2  =  1 

1 .0000 

-.5407 

-.4136 

0.5908 

0.0107 

-.0274 

FILTER 

0RDER--N1=2,N2=2 

1  .0000 

-.5572       0.0111 

-.4125 

0.6135       -.0404 

0.0099 

-.0327       0.0130 

FILTER 

0RDER--N1=2,N2=3 

1 .0000 

-.5574       0.0070 

0.0307 

-.4349 

0.6126       0.0007 

-.0287 

0.0284 

-.0398       -.0302 

0.0699 

FILTER 

ORDER— N1=3,N2  =  0 

1  .0000 

-.2071 

0.0564 

-.0402 

FILTER 

ORDER  —  N1=3,N2=1 

1.0000 

-.5402 

-.4132 

0.5903 

0.0074 

-.0216 

-.0043 

-.0133 

FILTER 

ORDER  —  N1  =  3,N2  =  2 

1.0000 

-.5571       0.0139 

-.4126 

0.6129       -.0406 

0.0073 

-.0206       0.0082 

-.0052 

-.0188       0.0113 

FILTER 

ORDER  — N1  =  3,N2  =  3 

1.0000 

-.5570       0.0082 

0.0302 

-.4365 

0.6125       0.0012 

-.0285 

0.0283 

-.0300       -.0469 

0.0781 

-.0130 

-.0126       0.0227 

-.0190 

Figure   14.      System  models  for  Eq.136  with  a  white  noise  sequence  input 
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or 


S(z1,22)=(F[r(«l,«2)])2  (137) 


S{zuz2)  =  (FlB(n\,n2)  *  e(nl,n2)lf  (138) 


or 


S(z„z2)  =  (B(z„  z2)E(zx,  z2))2  (139) 

If  the  input  is  a  unit  variance  white  noise  sequence,  the  power  spectrum  be- 
comes 

S{zuz2)  =  B\z„z2)    .  (140) 

We  can  now  use  the  all-pole  system  model  to  arrive  at  an  estimate  of  the  power  spec- 
trum 

S{zl,z2)  =  (       }        V- X- (141) 

where  A(zl,z2)  is  the  all-zero  prediction  error  filter  previously  derived. 

We  have  considered  three  examples  of  two-dimensional  all-pole  power  spectrum 
estimation.  In  all  examples  the  data  set  size  used  to  compute  the  autocorrelation  func- 
tion is  32  x  32.  In  the  first  example  a  sequence  consisting  of  two  cosines  at  normalized 
frequencies  (0.5,0.25)  and  (0.75,0.75)  is  corrupted  by  two-dimensional  unit  variance 
white  noise.  Figure  16  shows  the  surface  plot  of  the  estimated  power  spectrum  with  a 
4x4  filter  mask.  The  ideal  response  would  be  two  impulse  functions  at  the  appropriate 
frequencies.  The  plot  of  Figure  16  shows  several  spurious  peaks  along  a  ridge  which  is 
illustrated  by  the  contour  plot  in  Figure  17.  It  can  clearly  be  seen  from  the  surface  and 
contour  plots  in  Figures  18  and  19  that  the  estimation  of  the  power  spectrum  improves 
rapidly  with  an  increased  filter  size.  Figures  20  and  21  show  the  results  using  an  8  x  8 
filter.  The  surface  plot  in  Figure  20  clearly  shows  two  discrete  peaks  in  the  power 
spectrum.  Notice  that  there  is  a  slight  bias  in  the  estimates. 

In  the  second  example  three  cosines  at  frequencies  (0.5,0.25),  (0.75,0.75),  and 
(0.25,0.75)  are  used  to  generate  the  data.  The  surface  and  contour  plots  of  spectral  es- 
timates are  presented  in  Figures  22  through  27.  The  filter  masks  are  square  with  sizes 
4  x  4,  6  x  6,  and  8x8.    The  ability  of  these  filters  to  distiguish  the  spectral  peaks  again 
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Figure  15.      Generating  model  for  spectral  estimates 

improves  as  the  filter  size  increases  from  4  x  4  to  8  x  8.  While  the  4  x  4  filter,  as  shown 
in  Figure  22,  provides  very  little  useful  information,  the  8  x  8  filter  in  Figure  26  clearly 
indicates  three  distinct  peaks. 

The  third  example  investigates  the  capability  of  the  algorithm  to  differentiate 
between  two  closely  spaced  sinusoids  at  frequencies  (0.4375,0.4375)  and  (0.5625,0.5625). 
As  observed  in  Figures  28  through  31  the  algorithm  has  not  been  able  to  distinguish 
between  the  two  peaks.  However,  with  a  filter  size  8x8  (see  Figures  32  and  33),  the 
spectral  peaks  become  distinct  even  though  there  are  some  weak  spurious  peaks  present. 

The  examples  which  are  presented  in  Figures  16  to  33  show  the  usefulness  of  the 
algorithm  in  estimation  of  the  power  spectrum.  The  ability  of  the  algorithm  to  identify 
discrete  peaks  in  the  power  density  spectrum  increases  rapidly  with  filter  order.  As  can 
be  seen  in  Figures  32  and  33,  even  closely  spaced  peaks  are  readily  identified  by  an 
eighth-order  filter.  This  filter  requires  the  computation  of  64  coefficients  which  is  much 
simpler  than  using  Fourier  transforms.  All  plots  are  scaled  from  0  to  1  which  represents 
the  interval  0,7r. 

3.     Signal  Synthesis 

The  reflection  coefficients  for  a  single  frequency  data  set  were  determined  and 
used  to  implement  a  computer  simulation  of  the  synthesis  lattice.  The  resulting  lattice 
was  then  driven  by  a  pseudo  random  white  noise  sequence.  The  spectrum  of  the  output 
from  the  lattice  is  shown  in  Figure  34  for  a  4  x  4  structure.  Considering  the  small  size 
of  the  filter  and  the  poor  quality  of  the  input  noise  sequence,  the  results  are  acceptable. 
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4x4  FILTER 


u 


vo 


0.6 

Taxis'4     o.2 


oo 


0.2       °^^S 


Figure   16.      Surface  plot  of  spectral  estimate:     Y(nl,n2)=  cos(nl  n/2  +  n2  n;4)  + 
cos(nl  n  3/4  +  n2  7r  3/4) 
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4x4  FILTER 


i  r 


to 


-e- 


Figure   17.       Contour  plot  of  spectral  estimate:     Y(nl,n2)=  cos(nl  n',2   +   n.2  tt/4) 
+  cos(nl  n  3/4  +  n2  n  3/4) 
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6x6  FILTER 


o 


««« 


0   0 


Figure   18.      Surface  plot  of  spectral  estimate:     Y(nl,n2)=  cos(nl  n',2  +  n2  nlA)  4- 
cos(nl  n  3/4  4-  n2  n  3/4) 
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6x6  FILTER 


i 1 r 


"i 1 1 r 


to 
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0.2 


0.4  0.6 

$1   AXIS 


O.B  1.0 


Figure   19.       Contour  plot  of  spectral  estimate:     Y(nl,n2)=  cos(nl  7r/2   +   n2  tt/4) 
+  cos(nl  n  3/4  +  n2  n  3/4) 
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8x8  FILTER 


0  0 


Figure  20.      Surface  plot  of  spectral  estimate:     Y(nl,n2)=  cos(nl  rc/2  +  n2  n'4)  + 
cos(nl  n  3/4  +  n2  n  3/4) 
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8x8  FILTER 
i 1 1 1 1 r 


i r 


CO 


-e- 


J L 


J L 


0.2 


0.4  0.6 

4>1   AXIS 


0.8 


1.0 


Figure  21.       Contour  plot  of  spectral  estimate:     Y(nl,n2)  =  cos(nl  n/2  +   n2  tt/4) 
+  cos(nl  n  3/4  +  n2  n  3/4) 
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4x4  FILTER 


3 

2 


.    'ofr'i** 


o  0 


Figure  22.      Surface  plot  of  spectral  estimate:     Y(nl,n2)=cos(nl  re/2  +  n2  re/4)  + 
cos(nl  re  3,4  +  n2  re  3/4)  +  cos(nl  re/4,  n2  re  3/4) 
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Figure  23.       Contour  plot  of  spectral  estimate:     Y(nl,n2)=  cos(nl  n',2  +   n2  n,'A) 
+  cos(nl  n  3/4  +  n2  n  3/4)  +  cos(nl  tt/4,  n2  n  3/4) 


48 


6x6  FILTER 


o 


0  0 


Figure  24.      Surface  plot  of  spectral  estimate:     Y(nl,n2)=cos(nl  jt/2  +  n2  n;A)  + 
cos(nl  7T  3/4  +  n2  n  3/4)  +  cos(nl  rc/4,  n2  tt  3/4) 
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6x6  FILTER 


to 


-O- 


i  r 
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i  r 


<3^ 


J I  l I 


0.2 


0.4-  0.6 

♦  1   AXIS 


0.8 


1.0 


Figure  25.       Contour  plot  of  spectral  estimate:     Y(nl,n2)=cos(nl  nil  +   n2  n!A) 

+  cos(nl  n  3/4  +  n2  n  3/4)  +  cos(nl  njA,  n2  n  3/4) 
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8x8  FILTER 


120  - 


13 


*V 


,\  •** 


Figure  26.      Surface  plot  of  spectral  estimate:     Y(nl,n2)=  cos(nl  re/2  +  n2  re/4)  + 
cos(nl  re  3/4  +  n2  re  3/4)  +  cos(nl  re/4,  n2  re  3/4) 
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8x8  FILTER 
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Figure  27.       Contour  plot  of  spectral  estimate:     Y(nl,n2)  =  cos(nl  nil  +  n2  n,'4) 
+  cos(nl  n  3/4  +  n2  n  3/4)  +  cos(nl  tt/4,  n2  re  3/4) 
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4x4  FILTER 


1.0 


0  4  °'6 

0.2  Q'\y  Ax\S 


Figure  28.       Surface   plot   of  spectral   estimate:     Y(nl,n2)=cos(nl    rc7/16    +    n2 
77:7/16)  +  cos(nl7r9,16  +  n2  rr  9/16) 


53 


4x4  FILTER 


0.2 


0.4  0.6 

4>1   AXIS 


Figure  29.       Contour  plot  of  spectral  estimate:     Y(nl,n2)=  cos(nl    7r7/16    +    n2 
7i7,i6)  +  cos(ii1tt9;16  +  n2  n  9/16) 
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6x6  FILTER 


o 


*1  AX\S 


Figure  30.       Surface   plot   of  spectral   estimate:     Y(nl,n2)=cos(nl    nl/\6    +    n2 
7T7/16)  +  cos(nl7r9/16  +  n2  n  9/16) 
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6x6  FILTER 


0.4  0.6 

♦  1   AXIS 


Figure  31.       Contour  plot  of  spectral  estimate:     Y(nl,n2)  =  cos(nl    7t7/ 16    +    n2 
717/16)  +  cos(nl7r9/16  +  n2  tt  9/16) 
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8x8  FILTER 


< 


11  AXIS 


Figure  32.       Surface   plot   of  spectral   estimate:     Y(nl,n2)=cos(nl    ri7/16    +    n2 
tt7/16)  +  cos(nl7r9;i6  +  n2  ;r  9,  16) 
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8x8  FILTER 


0.4  0.6 

4>1  AXIS 


1.0 


Figure  33.       Contour  plot  of  spectral   estimate:     Y(nl,n2)  =  cos(nl    tt7,16    +    n2 
tt7,  16)  +  cos(nl7i9;16  +  n2  n  9/16) 
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SYNTHESIZED  SIGNAL 


o 


Figure  34.      Spectrum  of  signal  generated  from  white  noise:     Y(nl,n2)=cos(nl  tt/2 

+  n2  niA) 
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IV.     CONCLUSIONS 

The  linear  prediction  approach  to  analysis  of  stationary  one-dimensional  signals  has 
direct  applications  in  the  study  of  two-dimensional  signals.  The  solution  of  the  2-D 
normal  equations  in  terms  of  orthogonal  polynomials  provides  an  efficient  and  effective 
method  for  system  modeling  and  power  spectral  density  estimation.  There  also  exist 
two-dimensional  lattice  structures  that  are  analogous  to  their  one-dimensional  counter- 
parts. These  lattices  can  be  used  as  an  alternate  implementation  of  prediction  error  fil- 
ters and  signal  syntcsis  filters. 

The  simple  nature  of  the  recursive  solution  to  the  2-D  normal  equations  developed 
in  this  thesis  provides  an  excellent  method  for  obtaining  solutions  to  problems  that  are 
not  dependent  on  the  stability  of  the  derived  filter.  For  example,  system  modeling 
achieves  very  accurate  results  assuming  the  system  being  modeled  is  stable.  Signal  syn- 
thesis on  the  other  hand  can  not  be  considered  reliable  since  there  is  nothing  to  guar- 
antee a  stable  filter.  This  can  be  attributed  to  the  lack  of  autocorrelation  matching  in 
two-dimensional  signals.  An  excellent  discussion  of  this  property  can  be  found  in  [Rcf. 
4J.  Simply  stated,  there  arc  more  independent  autocorrelation  coefficients  than  there  are 
coefficients  in  the  filters  that  we  have  found.  We  have  ignored  the  second-  and  fourth- 
quadrant  coefficients  which  arc  independent  from  the  first-  and  second-quadrant 
cocfliccints  we  have  been  using. 

Considering  the  computational  expense  of  other  spectrum  estimation  techniques  the 
algorithm  presented  in  this  paper  is  valuable  if  used  only  for  that  purpose.  The  problem 
of  signal  synthesis,  while  more  difficult,  was  met  with  some  success  as  evidenced  by 
Figure  34.  This  limited  success  invites  further  study  into  autoregrcssive  modeling  of  2-D 
sienals. 
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